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We consider the Kuramoto model of an ensemble of interacting oscillators allowing for an arbitrary 
distribution of frequencies and coupling strengths. We define a family of traveling wave states as 
stationary in a rotating frame, and derive general equations for their parameters. We suggest 
empirical stability conditions which, for the case of incoherence, become exact. In addition to 
making new theoretical predictions, we show that many earlier results follow naturally from our 
general framework. The results are applicable in scientific contexts ranging from physics to biology. 



Almost every real-life physical system involves a large 
number of interacting subsystems. In many cases, they 
can be treated as a population of interacting phase os- 
cillators that can be described in terms of the Kuramoto 
model (KM) which we consider in the form 
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Here 0i, u>i and Ki are respectively the ith oscillator's 
phase, natural frequency, and strength of coupling to the 
other oscillators; and uji and Ki are randomly chosen 
from a probability density g(ui,K). The KM has been 
used in a variety of applications, ranging from brain dy- 
namics and human crowd behavior to Josephson junction 
arrays and neutrino flavor oscillations 0t3; so th & t the 
analysis of its dynamics is of high topical interest and 
broad applicability in science. 

Hitherto, KM solutions have been obtained only for 
particular cases of g(uj, K) . For the classical KM I;] , Ki is 
fixed, so g(u>,K) — g(u))8(K — Kq), with a restriction for 
g(u>) to be unimodal and symmetric about its maximum. 
Other works ((J 0] go beyond the classical KM by consid- 
ering other particular forms of g(u,K), but no attempt 
has been made to solve ([1} for arbitrary g(ui,K). Here 
we take a major step towards filling this gap, deriving 
general equations and approximate stability conditions 
for the family of traveling wave (TW) states, defined as 
those that are stationary in some rotating frame. 

The collective behavior of Kuramoto oscillators is de- 
scribed by the order parameter 



Z = Re^ = — V 
N ^ 



(2) 



The mean field strength R quantifies the synchronisation 
between the oscillations and is usually the main quantity 
of interest. Substitution of ([2]) into (pj yields 



= u>i — KiR sin(#i — ip). 



(3) 



In the continuum limit, N — > oo, (O is treated using 
the probability density function (PDF) f(9,u>,K,t), i.e. 
the probability that an oscillator has phase 6, coupling 
strength K and frequency uj at time t. In many cases [8[ 



the PDF may be represented using the OA ansatz Q 
f(0, u , K, t) = ?^p- J 1 + J2 K, t)] n 
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where c.c. denotes complex conjugate, and a satisfies 
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Z = Re 1 ^ 



a*{uj,K,t)g{uj,K)dujdK. (6) 



The integrals are taken over (—00,00) if unspecified. De- 
noting a(uj,K,t) = p(uj,K,t)e i 'l > ^' K ^, © becomes 
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(1 + p 2 ) sin(0 + -0) + up = -p 



df 



We now define some terms. The KM equations ([3]) do 
not change form when transformed to a frame rotating 
at £1 (9i — > 6i~VLt): this is equivalent to changing the fre- 
quency distribution g(u>, K) — > g(uj+fl, K) with uj always 
denoting the natural frequency in the current rotating 
frame. Thus, all rotating frames are physically equiva- 
lent. By a stationary state (SS) we will mean a state with 
a time-independent PDF: df /dt = 0. For an SS, it is eas- 
ily shown that Z remains constant. However, the SS is 
stationary only in its own rotating frame (except in the 
case of incoherence) because it becomes time-dependent 



when changing frame (Z 



-nt\ 



Thus, apart from 



its order parameter and stability, an SS is also charac- 
terized by its frame frequency, i.e. the frequency ft of the 
rotating frame in which it is stationary. 

The frame with zero mean frequency (uj) 
I J u)g(uj,K)diu/,dK = will be called the natural 
frame, and we will understand g(uj,K) as the distribu- 
tion in the natural frame; Q will denote the current frame 
frequency with respect to natural one. We will see that 
the usual TW states are in fact stationary, but just not 
in the natural frame. So we introduce the notion of a 
family of TW states that are stationary in some rotat- 
ing frame and which can be described by the OA ansatz. 
This family also includes states stationary in the natural 
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frame (Q = 0), such as partially synchronized states and 
the 7r-state. We will call these kinds of TW states natural 
states. For convenience we introduce the notation 



9± = g(±u + ^, K), g = g{uj, K). 



We call a distribution uncorrelated if the distributions of 
u> and K are independent (so g(uj,K) = g(u))T(K)) and 
correlated otherwise. We call it symmetric if g(w, K) — 
g(—uj,K) and asymmetric otherwise. 

For stationary states §f = ^ = 0. Substituting this 
into ([7]) and noting that p < 1 for the OA ansatz to be 
valid, one can show that all stable SSs arc of form 



a(w, K)e 



■ u)—sign(u))^/uj 2 — K 2 R 2 
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if U > \K\R 



if U < \K\R. 
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There can be a minus before the V ' K 2 R 2 — to 2 for \u\ < 
\KR\, but it always corresponds to an unstable position 
on the phase circle. We are interested only in stationary 
states that can be realized, so we discard this solution. 
Using (8) in one can reconstruct the full PDF 

f(0,ui,K) = g(u+Sl,K)x 

j 5{6-iP- arcsin( + ttR{-K)) if \u\ < \K\R 

{ Vu 2 -K 2 R 2 /2tt • r I I. i MD 
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where H(-) is a Heaviside function. Eq. (fTUf is well 
known for constant K > [10], while for K < it cor- 
rectly reflects a change of the stable phase difference to 
it. Note that, for a bimodal delta-distribution of K, the 
difference between the complex phases of oscillators with 
K > and K < (arg Y^i-K^o e%di ) ec l ua l s n omv f° r 
the natural state, and differs from it in the TW state 
with fl ^ 0- This might appear inconsistent with 
(fTOj) . but it is not: (| 10[) implies a difference of it between 
the complex phases of populations with K ^ only if 
g(oj + H,K) = g{—uj + £l,K) which cannot be true for 
f2 7^ 0. The representations and (jTOJ) are equivalent. 
Each of the following expressions can be derived from ei- 
ther, but use of the OA ansatz stationary solution (J9j) 
enormously simplifies derivations. 

Using R = J a*(uj,K,t)e- i,f 'g{uj,K)dujdK. Tak- 
ing its real and imaginary parts and using yields: 

+ \K\R 

F R (R, n)^-R + J^ J gWK 2 R 2 -u; 2 dio = 0, 
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F n (R,n)= I I Lug+dcu- 



J [g+ - g-] V^ 2 - K 2 R 2 du} = 0. 
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These are the general self-consistency conditions (SCC), 
for any g(u), K). They determine the mean field strength 



and frame frequency of the SS. Note that for symmetric 
frequency distributions, g(uj,K) = g(—uj,K), one solu- 
tion of the second equation is always O = 0, and that 
TW states are born in pairs of similar R and opposite 
n f (|11[) is invariant under f2 — > — O). For asymmetric 
frequency distributions, the second of (jlip implies indi- 
vidual TW states, not pairs. 

Having determined the parameters of the stationary 
state, we need to find the conditions for its stability. In 
the general case, this is an extremely challenging prob- 
lem. To tackle it, we can derive an approximation that 
we will call the empirical stability condition (ESC). We 
postulate that the stability of a particular SS should be 
connected with the stability of its SCC (TTT|) . which may 
be written as Fr + iFn = 0. 

Next, we consider the stationary state to be stable if 
its SCC is stable with respect to R in the real direction 
and £1 in the imaginary one, i.e. stability is equivalent to 
< 0, ^ < 0. Using (jnj, our ESC are then 



\K\R 




ui 2 g+dui 



-\k\r K 2 R 2 \/K 2 R 2 - , 

u[g++g-\ duj 
\k\r Vuj 2 - K 2 R 2 



< 1, 



(12) 



< 0. 



The first of (TT3]) characterizes stability with respect to 
changes in mean field amplitude, while second refers to 
creation of a TW, i.e. to a change of frame. 

Using Eqs. (1121) . one can obtain a universal condition 
for the stability of incoherence. The 1st of ([12")) for R — » 
reduces to ^ J Kg(fl, K)dK < 1. Note however that, for 
R = 0, Eqs. (fTTj) arc satisfied for any ft (and the 2nd of 
Eqs. p^|) is zero). This reflects the physically self-evident 
fact that an incoherent state must be the same in any 
rotating frame. Thus, the main question is now what f2 
to use for the determination of stability. It seems most 
appropriate to use the O of a frame in which clustering 
can evolve from incoherence, i.e. in which Fq(R, Q) — 
is satisfied up to next order of R as R — > 0. For a 
distribution continuous over lj, our ESC for incoherence 
becomes 



KdK 



g(oj + fiW , K) -g(-w + fiW , K) 



duj = 0, 



: J Kg(fl {t) ,K)dK < ^ 



(13) 

where in the inequality we choose the maximum over the 
solutions of the first of p3]) : it is clear that, for 
overall stability, incoherence should be stable in all its 
frames. We now observe that (| 131) is in fact the exact 
stability condition, or at least the transition point. 

This result can be proven more rigorously. A perturba- 
tion to the stationary state should grow according to (|5|) . 
Using the standard procedure [111 ] with the OA ansatz in- 
stead of the full PDF, we linearize ([5]) about incoherence 
(for which a = 0), substituting a(ui, K) = + j3(u>, K)e xt 
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to yield (A + ioj)/3(oj, K) = (K/2) J (3(w,K)g + dujdK. 
Denoting B = J /3(uj, K)g + dujdK, we obtain the self- 
consistent form f3(u,K) = -jJTIu' Substituting the lat- 
ter into the definition of B itself (note that K then en- 
ters the integral), we obtain the self-consistency equation 
J ^xTfaj— — 2, from which the perturbation growth ex- 
ponent A can be determined. The state becomes linearly 
unstable when ReX exceeds zero. To find the transition 
point we substitute A = Ax + 1X2 into the self-consistency 
equation and take the limit Ai — > + , yielding 



1 = | J Kg(uj + n- \ 2 ,K)dK, 
Kg(u + Q - \ 2 ,K)duidK 



= 



(14) 



This is completely equivalent to the transition point pre- 
dicted by (IT3"1) . The only difference is that here the fiW 
are determined through A 2 . This justifies the intuition 
underlying the ESC. However, (fl"4"|) determines the tran- 
sition point when incoherence becomes linearly unstable, 
and it does not predict the inequality in (fT3"]) . It might 
also fail in rare cases when the OA ansatz is invalid. 

In the simplest case of an uncorrelated g(ui, K) = 
g{uj)Y(K) with g(co) continuous, unimodal and symmet- 
ric about its single peak, one has fiW = in 03]). Thus, 
we obtain for this very general case the simple result 



(K) c 



(15) 



where (K) c is the critical value where incoherence be- 
comes unstable. Eq. (fT5)) reduces to the Kuramoto result 
K, 



in the case of constant K, and to the recent 
p} (Eq. 12) for the case studied there. 



results of 

In the case of a bimodal Lorenzian frequency distribu- 



tion g(ui) 



2tt 1 



(oj—u 



(p]l yields 



(K) c = 



2(w 2 + 7 2 )/ 7 if w < 7 
47 if w > 7 



(16) 



which includes corresponding results of [6]. For uncor- 
related distributions, Eqs. (|13|) imply that, for any dis- 
tribution of couplings T(K), incoherence stability is de- 
termined only by the mean value (K) and not by higher 
moments, being very similar to the constant K case. 

Another interesting case involves correlated distribu- 
tions. For example, even if the conditional distribu- 
tion of frequencies g(u>\K) is unimodal and symmet- 
ric, there might be solutions of (fT3| other than 
zero, because the corresponding integral might vanish af- 
ter integration over K. Consider e.g. the distribution 
g(u,K) = [(1 - p)S(K K X ) + P S(K - K 2 )] 7J ^} E) 

Denoting 71,2 = 7(^1,2) and using (IT^|) . one can ob- 
tain an elegant condition for the stability of incoherence 
max |(i^i + £^ ; J|L if ( K )(K/j 2 (K))<0} < 2. 

Thus a new maximum will appear only when K\ 2 has 
both signs. Note that, for continuous uncorrelated dis- 
tributions, incoherence is always stable if (K) < 0; but 
for some kinds of correlated distributions it can be unsta- 
ble even in this case, giving rise to synchronized states. 

In certain cases the above formulae can be reduced to 
a simpler form. Consider a multimodal Lorenzian fre- 
quency distribution in its most general form 



g(u,K)=T(K)J2 



\ (W - W„) 2 +7n 



Qn(K) = 1VK 

71=1 

where, additionally, the parameters are allowed to de- 
pend on K: q n (K), j n (K), uj n (K). In this case we inte- 
grate (fTTj) over to explicitly [l2[ and obtain 
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n=l J 



KR \ hl ^2 
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q n T(K) 
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dK l n 



V2 



^Jy/{K*B? + 7 2 - + 4 7 2 ^ + {K 2 R 2 + 7 2 -02)} 

^/{K*R? + H - + 4 7 |fi2 _ (K 2 R 2 + 7 2 _ Q2 )} 



= 

= 
(18) 



where O n = Q— u) n (K ). These expressions for determina- 
tion of R and £1 become fully algebraic when one has mul- 
timodal delta distributions of K . Thus, we can readily 
reproduce the results of @, 0] for R in the natural state. 
Also, (|18|) can be expanded over 7n to obtain approxi- 
mate expressions for more complicated T{K). Empirical 
stability conditions (fT2"j) may be calculated from (fTS]) as 



T§f - °' W - °- For incoherence stability, the first of 



([12]) takes the simple form J^n=i / 



Kr{K)g n n n _ 



0. 



Consider the case of a unimodal frequency distribution: 



g( i o,K) = T(K)-£S^ 



It immediately follows from 
fl8l) that, for existence of the TW state (solutions O ^ 0), 
one needs K of both signs. Eqs. (fT8|) and (fT2j) reduce to 



R 



< 



KR 



dK 



:{-i(K) + v/X 2 i? 2 + 7 2 (X)} 
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< Rf 



(19) 
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FIG. 1. (Color online) Mean field strength R and frame frequency Q for different distributions g(ui,K) as functions of 
parameter p oc (K) which characterizes the distribution of couplings: F2{K,p) = (1 — p)S(K — Ki) + pS(K — K-i) or 



-f-5(K - Ki) + ^5{K - K 2 ) + pS(K - K 3 ) (we also denote L(uj, 7) 



(a) g(u,K) = L(oj,~f)r 3 (K,p) 



with 7 = 0.025, = -1.5, K 2 = -0.75, K 3 = 2; (b) g(u,K) = L(u,i\K\)T 2 (K,p) with 7 = 0.1, Ki. = -0.5, #2 = 1; 
(c) g(u,K) = T 2 {K,p)e-" 2 l 2a2 /-^o with o = 0.05, /fi = -0.5, K 2 = 1; (d) g(u,K) = T 3 (K,p)\l + erf(a(w - 

Mo)/v / 2o-)]e" ( "" Mo)2/2CT2 /V^Ktr, where fj, = -V2va/y/w(l + a 2 ) (ensures (w) = 0) and a = 0.1, a = 2,^1 = -0.5, ^2 = 1. 
Theoretical predictions (lines) are compared with numerical simulations (circles). I = incoherence; NS = natural state; TW = 
traveling wave. Full lines imply stable states, according to ESC; dashed lines imply unstable. Short black lines indicate TW 
states for which ESC fails. In the insets of (a) and (b) TW states with Q, < are omitted for clarity, since they are the same, 
whereas in (d) both are shown. The inset of (c) shows the hysteretic region in more detail. The simulations used N = 25600 
oscillators and a 6th order Runge-Kutta algorithm with time step 0.01 s for 500 s. R and Q are averages over the last 100 s. 



which are fully algebraic in the case of T(K) = 
Y^=iPm5{K - K m ),J2 m P™ = 1 - For c onstant 
7 and K we get the well known R = \J 1 — ^ , 
while for the bimodal distribution (N p = 2) para- 
metric R 2 = pisign^xX^ i? 2 + (7/ifi) 2 - j/\Ki\) + 

P2Sign(ir 2 )(V R 2 + {l/K 2 ) 2 - -i/\K 2 \) 7 which may be 
brought to the form of Eq. (13) of [7]. 

The predictions of (|T3|) . (fT8)l . (fT9t are compared with 
numerical simulations in Figs. l(a,b) where a novel tri- 
modal distribution of K and correlated joint distribution 
are studied. In both cases the agreement is excellent. 
Note the interesting decrease in R with increasing p (and 
so (K)) in a region around p = 0.62 of Fig. 1(a). The 
validity of the general Eqs. (TTTT) - (p~5|) is tested in Fig. 1(c) 
for Gaussian and in (d) for asymmetric frequency distri- 
butions. From (c) it is clear that the theory can predict 
phenomena like the co-existence of states, and hysteresis, 
which can be hard to catch numerically. In Fig. 1(d) O 
seems to go to zero at the higher p but it does not quite 
get there, reflecting that there are only TW states. We 



also see that, according to ESC, there should exist one 
more stable TW state in both Fig. 1(c) and (d). Here 
ESC fails because, in practice, such states are unstable. 
Although ESC is not infallible, its good performance in 
so many cases suggests that some general stability condi- 
tions may be derived, and that an analysis of when ESC 
works and when it does not might be fruitful. 

Note that, except for (TT3"|) . the formulas were derived for 
the TW family of states that cannot describe purely non- 
stationary states like standing waves. There may perhaps 
exist other special states not included in our framework, 
e.g. some SSs not described by the OA ansatz. Never- 
theless, if such states exist it would be of considerable 
interest to find them, and our framework may greatly 
help through comparison with numerical simulations. 

In conclusion, our new framework for analysis of the 
KM with arbitrary g(uj, K) incorporates many previous 
results and enables us to predict new phenomena. Our 
findings may be useful in the diverse applications of the 
KM, e.g. those mentioned in the introduction, especially 
when the couplings are expected to be heterogeneous. 
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